## Measuring Networks in the Field
## Larson and Lewis
## Replication Code
## Updated 12-14-2018

#####################################################
## This code uses data from 8 measured networks.  Six are from:

##Banerjee, Abhijit;Chandrasekhar, Arun G.;Duflo, Esther;Jackson, Matthew O., 2013, "The Diffusion of Microfinance", https://hdl.handle.net/1902.1/21538, Harvard Dataverse, V9

## Two are from:

##Larson, Jennifer;Lewis, Janet, 2016, "Replication Data for: Ethnic Networks", https://doi.org/10.7910/DVN/W1TDGZ, Harvard Dataverse, V1


## This code simulates sampling from each of these eight networks, stores network statistics computed on the samples, and compares the rank of the statistics for the sample to the true network.  


library(igraph)
set.seed(12345)

## Import data, first 6 villages from Banerjee, Chandrasekhar, Duflo, and Jackson'ss microfinance data:

data1 <- read.csv("nwdata/adj_allVillageRelationships_vilno_1.csv", header=FALSE)
data2 <- read.csv("nwdata/adj_allVillageRelationships_vilno_2.csv", header=FALSE)
data3 <- read.csv("nwdata/adj_allVillageRelationships_vilno_3.csv", header=FALSE)
data4 <- read.csv("nwdata/adj_allVillageRelationships_vilno_4.csv", header=FALSE)
data5 <- read.csv("nwdata/adj_allVillageRelationships_vilno_5.csv", header=FALSE)
data6 <- read.csv("nwdata/adj_allVillageRelationships_vilno_6.csv", header=FALSE)

## And the 2 Ugandan villages from Larson and Lewis:

dataac <- read.csv("nwdata/ac_adj.csv", header=TRUE)
datamc <- read.csv("nwdata/mc_adj.csv", header=TRUE)

## Convert all eight to network objects

nw1 <- graph.adjacency(as.matrix(data1))
nw2 <- graph.adjacency(as.matrix(data2))
nw3 <- graph.adjacency(as.matrix(data3))
nw4 <- graph.adjacency(as.matrix(data4))
nw5 <- graph.adjacency(as.matrix(data5))
nw6 <- graph.adjacency(as.matrix(data6))
nwac <- graph.adjacency(as.matrix(dataac))
nwmc <- graph.adjacency(as.matrix(datamc))

rm(data1, data2, data3, data4, data5, data6, dataac, datamc)

######################################################
## Sample size, node-level attributes
######################################################


# Features of interest
nw_features <- c( "outdeg", #1
                 "neighb1", #2
                 "neighb2", #3
                 "neighb2only", #4
                 "eigen", #5
                 "transitivity" #6
                )


calculate_nw_features <- function(nw) {
  outdeg <- degree(nw, v= V(nw), mode = "out")
  neighb1 <- ego_size(nw, order = 1, nodes= V(nw), mode = "all")-1
  neighb2 <- ego_size(nw, order = 2, nodes= V(nw), mode = "all")-1
  neighb2only <- neighb2 - neighb1 
  eigen <- eigen_centrality(nw)$vector
  transitivity <- transitivity(nw, v = V(nw), type = "local")
  
  values <- as.matrix(rbind(outdeg, neighb1, neighb2, neighb2only, eigen, transitivity)) # add name of any new features here so that the function will return it
  
  return(values)
}

## Helper function: RMSD, assume length(a) = length(b)
rmsd <- function(a,b) {
  (sum((a-b)^2, na.rm=TRUE)/ length(a) )^(1/2)
}


### Main function to calculate cor and rmsd
# ranksims() constructs the output matrices, calculates network feature values, and calculates the Kendall correlation and RMSDs
# Calls calculate_nw_features() on the true, closed, and open networks
# Returns a list of 5 matrices: true/closed correlation, true/open correlation, true/closed RMSD, true/open TMSD, and the true values
# @param nw = Network object
# @param nsims = Num sims
# @param psamp = proportion to sample for open/closed
# @param feature_names = Character vector of names of network features of interest (used for labeling matrices etc) -- will most likely pass nw_features to this argument
# @param seed = seed

ranksims <- function(nw, nsims, psamp, feature_names, seed = 12345) {
  set.seed(seed)
  
  adj <- as.matrix(get.adjacency(nw))
  
  # Batch instantiate matrices to hold correlation and RMSD vals for closed/open, respectively
  nrows <- length(feature_names)
  
  for (matrix_name in c("ccor", "ocor", "crmsd", "ormsd")) {
    assign(matrix_name, matrix(NA, nrow = nrows, ncol = nsims))
  }
  
  rownames(ccor) <- feature_names
  rownames(ocor) <- feature_names
  rownames(crmsd) <- feature_names
  rownames(ormsd) <- feature_names
  
  # Calculate the nw features for the full/true network
  true_nw_features <- calculate_nw_features(nw)
  
  # Sample the network
  nkeep <- ceiling(psamp * vcount(nw))
  
  for(i in 1:nsims) {
    samp <- sample(1:vcount(nw), size = nkeep, replace = FALSE)
    
    # Sample closed network
    tempadjc <- adj[samp, samp]
    tempc <- graph.adjacency(tempadjc)
    
    # Sample open network
    tempadjo <- adj
    totoss <- setdiff(1:nrow(adj), samp) #find the unsampled rows 
    tempadjo[totoss,] <- 0 #toss them (keeping unsampled columns)
    tempo <- graph.adjacency(tempadjo)
    
    ## Calculate true rank for the sampled vertices
    true_sampled <- true_nw_features[, samp]  # Get subset of the full matrix of true values: True values for sampled vertices # Note: This is Features x Sampled Vertices
    # Calculate rank
    true_ranks_sampled_Vs <- t(sapply(nw_features, function(x) rank(true_sampled[x, ]))) # Note: transposed matrix to maintain Features x Vertices format
    
    ## Calculate network features, and then rank, for the closed/open samples
    # Closed:
    closed_nw_feature_vals <- calculate_nw_features(tempc)
    closed_ranks <- t(sapply(nw_features, function(x) rank(closed_nw_feature_vals[x, ]))) # transposed here too
    
    # Open:
    open_nw_feature_vals <- calculate_nw_features(tempo)
    open_nw_feature_vals <- open_nw_feature_vals[, samp]
    open_ranks <- t(sapply(nw_features, function(x) rank(open_nw_feature_vals[x, ]))) # transposed here
    
    ## Calculate and store correlations and RMSD
    for (n in 1:nrows) { 
      current_feature <- nw_features[n] # nrows = num of network features we're interested in. Iterate through the network features, using this to grab the column of rankings from the true rank matrix and the closed or open rank matrix
      
      ## Calculate corr: 
      # Select the true, open, and closed rankings
      current_true_ranks <- true_ranks_sampled_Vs[current_feature, ]
      current_closed_ranks <- closed_ranks[current_feature, ]
      current_open_ranks <- open_ranks[current_feature, ]
      
      # Calculate Kendall rank cor for closed/true and open/true, store in ccor and ocor matrices
      ccor[n, i] <- cor(current_true_ranks, current_closed_ranks, method = "kendall") # Reminder: i = simulation num, n = feature num
      ocor[n, i] <- cor(current_true_ranks, current_open_ranks, method = "kendall")
      
      ## Calculate rmsd: 
      #Select the true, open, and closed values
      current_true_vals <- true_sampled[current_feature, ]
      current_closed_vals <- closed_nw_feature_vals[current_feature, ]
      current_open_vals <- open_nw_feature_vals[current_feature, ]
      
      # Calculate RMSD for closed/true and open/true, store in crmsd and ormsd matrices
      crmsd[n, i] <- rmsd(current_true_vals, current_closed_vals)
      ormsd[n, i] <- rmsd(current_true_vals, current_open_vals)
    }
    
    if (i%%50 == 0) {
      print(i) # sim counter    
    }
    
  }
  
  # Create list of matrices to return
  out <- list("ccor" = ccor, "ocor" = ocor, "crmsd" = crmsd, "ormsd" = ormsd, "truth"=true_sampled)
  return(out)
}



# Example output with a single sample size
nw1_sims <- ranksims(nw1, 1000, 0.2, nw_features)

# This generates a list of 5 matrices
# nw1_sims$ccor  Get matrix of Kendall rank correlation matrix for values from the true and closed sampled network 
# nw1_sims$ocor  Get matrix of Kendall rank correlation matrix for values from the true and the open sampled network
# nw1_sims$crmsd  Get matrix of RMSD for values from the true and closed sampled network
# nw1_sims$ormsd  Get matrid of RMSD for values from the true and open sampled network
# nw1_sims$truth  Get matrix of true feature values for the sampled nodes. 

# Note: In all matrices, rows correspond to the network features named in nw_features 
# In $truth matrix, columns are labeled with the original node number reflecting which nodes were selected when sampled. In all of the other matrices, columns are labeled with 1:num_sampled_nodes


##############################################################################
### Now run sims for multiple sample sizes

## Helper function to generate rank_sims for each value in a list of p_sample_sizes
# Takes all of the @params for ranksims() and passes them to ranksims(), except that it also 
# takes a vector of p_sample_sizes and asks for a descriptive string (nw_name) to give to the return value
# @param nw = Network object
# @param nw_name = String name to give to the return values
# @param nsims = Num sims
# @param ps = Vector of p sample sizes 
# @param feature_names = Character vector of names of network features of interest
# @param seed = seed

generate_ranksims <- function(nw, nw_name, nsims, ps, feature_names, seed = 12345) {
  sim_list <- c() # list of list of matrices to return
  
  for (p in ps) {
    sim_name <- paste(nw_name, p, sep = "_") # this will construct a name for the specific simulation based on the nw_name and the p_size
    assign(sim_name, ranksims(nw, nsims, p, feature_names)) # create sim for this p_size
    sim_list <- c(sim_list, list(get(sim_name))) # add ranksims() result to the list
  }

  return(sim_list)
}

# Set values to use throughout
nsims <- 1000
ps <- c(.1, .2, .3, .4, .5, .6, .7, .8, .9)

nw1_sims_list <- generate_ranksims(nw1, "nw1", nsims, ps, nw_features, seed = 12345)
nw2_sims_list <- generate_ranksims(nw2, "nw2", nsims, ps, nw_features, seed = 12345)
nw3_sims_list <- generate_ranksims(nw3, "nw3", nsims, ps, nw_features, seed = 12345)
nw4_sims_list <- generate_ranksims(nw4, "nw4", nsims, ps, nw_features, seed = 12345)
nw5_sims_list <- generate_ranksims(nw5, "nw5", nsims, ps, nw_features, seed = 12345)
nw6_sims_list <- generate_ranksims(nw6, "nw6", nsims, ps, nw_features, seed = 12345)
nwac_sims_list <- generate_ranksims(nwac, "nwac", nsims, ps, nw_features, seed = 12345)
nwmc_sims_list <- generate_ranksims(nwmc, "nwmc", nsims, ps, nw_features, seed = 12345)



#######################################################
## Plotting the results of sampling simulations



## Helper function to plot the Kendall rank correlation for a given network and a given network feature
# @param sim_list_name : Name of the list holding the list of simulations for p_sample_sizes
# @param nw_name : Name of network as a string for labelling
# @param feature_name Name of feature as a string (choose from options in nw_features )
# @param nsims : nsims from before
# @param ps : ps from before
plot_feature_cor <- function(sim_list, nw_name, feature_name, nsims, ps) {
  title <- paste(as.character(feature_name), "in", as.character(nw_name), sep = " ")
  plot(x = rep(.1, nsims), y = sim_list[[1]]$ccor[feature_name,], xlim = c(0,1), ylim = c(0,1), xlab = "Sample Size", ylab = "Kendall Rank Correlation", main = title)
  
  points(x = rep(.11, nsims), y = sim_list[[1]]$ocor[feature_name,], col = "skyblue")
  for(i in 1:length(ps)) {
    points(x = rep(ps[i], nsims), y =sim_list[[i]]$ccor[feature_name,])
    points(x = rep(ps[i]+.01, nsims), y = sim_list[[i]]$ocor[feature_name,], col = "skyblue")
  }
}



## Function to generate a plot of Kendall rank cor for all network features for a given network
# Calls plot_feature_cor() on each of the network features for the given network
# @param sim_list Name of list of simulations for p_sample_sizes
# @param nw_name String name to give this network
# @param nw_features List of names of network features of interest
# @param nsims Nsims
# @param ps Vector of p_sample_sizes
#
plot_cor_all_features <- function(sim_list, nw_name, nw_features, nsims, ps) {
  par(mfrow=c(2,3), oma=c(0,0,2,0)) # change this to fit more plots
  
  for (i in 1:length(nw_features)) {
    plot_feature_cor(sim_list, nw_name, nw_features[i], nsims, ps)
  }
  title(paste("Kendall Rank Correlations for ", nw_name, sep = " "), outer=TRUE)
  par(mfrow=c(1,1)) # reset plot params
}

## Plots all six features per network 

## SI Figure 1:
pdf("SIf1.pdf")
plot_cor_all_features(nw1_sims_list, "NW1", nw_features, nsims, ps)
dev.off()

## SI Figure 2:
pdf("SIf2.pdf")
plot_cor_all_features(nw2_sims_list, "NW2", nw_features, nsims, ps)
dev.off()

## SI Figure 3:
pdf("SIf3.pdf")
plot_cor_all_features(nw3_sims_list, "NW3", nw_features, nsims, ps)
dev.off()

## SI Figure 4:
pdf("SIf4.pdf")
plot_cor_all_features(nw4_sims_list, "NW4", nw_features, nsims, ps)
dev.off()

## SI Figure 5:
pdf("SIf5.pdf")
plot_cor_all_features(nw5_sims_list, "NW5", nw_features, nsims, ps)
dev.off()

## SI Figure 6:
pdf("SIf6.pdf")
plot_cor_all_features(nw6_sims_list, "NW6", nw_features, nsims, ps)
dev.off()

## SI Figure 7:
pdf("SIf7.pdf")
plot_cor_all_features(nwac_sims_list, "NW7", nw_features, nsims, ps)
dev.off()

## SI Figure 8:
pdf("SIf8.pdf")
plot_cor_all_features(nwmc_sims_list, "NW8", nw_features, nsims, ps)
dev.off()

## Now plot mean by sample size for open for each feature:

## First make vectors of mean correlations by sample size:

## Create a matrix for each network that stores mean correlation b/t feature rank in open sampling and true rank among those sampled, where rows are feature and columns are sample size

nw1_mat <- matrix(nrow = length(nw_features), ncol = length(ps), data = NA)
for (i in 1: length(ps)){
	nw1_mat[,i] <- apply(nw1_sims_list[[i]]$ocor, 1, mean, na.rm=TRUE)
}

nw2_mat <- matrix(nrow = length(nw_features), ncol = length(ps), data = NA)
for (i in 1: length(ps)){
	nw2_mat[,i] <- apply(nw2_sims_list[[i]]$ocor, 1, mean, na.rm=TRUE)
}

nw3_mat <- matrix(nrow = length(nw_features), ncol = length(ps), data = NA)
for (i in 1: length(ps)){
	nw3_mat[,i] <- apply(nw3_sims_list[[i]]$ocor, 1, mean, na.rm=TRUE)
}

nw4_mat <- matrix(nrow = length(nw_features), ncol = length(ps), data = NA)
for (i in 1: length(ps)){
	nw4_mat[,i] <- apply(nw4_sims_list[[i]]$ocor, 1, mean, na.rm=TRUE)
}

nw5_mat <- matrix(nrow = length(nw_features), ncol = length(ps), data = NA)
for (i in 1: length(ps)){
	nw5_mat[,i] <- apply(nw5_sims_list[[i]]$ocor, 1, mean, na.rm=TRUE)
}

nw6_mat <- matrix(nrow = length(nw_features), ncol = length(ps), data = NA)
for (i in 1: length(ps)){
	nw6_mat[,i] <- apply(nw6_sims_list[[i]]$ocor, 1, mean, na.rm=TRUE)
}

nw7_mat <- matrix(nrow = length(nw_features), ncol = length(ps), data = NA)
for (i in 1: length(ps)){
	nw7_mat[,i] <- apply(nwac_sims_list[[i]]$ocor, 1, mean, na.rm=TRUE)
}

nw8_mat <- matrix(nrow = length(nw_features), ncol = length(ps), data = NA)
for (i in 1: length(ps)){
	nw8_mat[,i] <- apply(nwmc_sims_list[[i]]$ocor, 1, mean, na.rm=TRUE)
}


## Now plot each of the 8 village's mean rank correlations between features and true rank
colors <- c("orchid", "violetred", "turquoise", "steelblue1", "palegreen", "sienna1")


plot(x = ps, y = nw1_mat[1,], col = colors[1], type = "l", ylim = c(0,1), xlab = "Sample Size", ylab = "Mean Kendall Rank Correlation", main = "NW1")
for (i in 2:length(nw_features)){
	lines(x = ps, y = nw1_mat[i,], col = colors[i])
	}
legend("bottomright", legend = nw_features, lty = c(1,1,1,1,1,1), col= colors)
	
plot(x = ps, y = nw2_mat[1,], col = colors[1], type = "l", ylim = c(0,1), xlab = "Sample Size", ylab = "Mean Kendall Rank Correlation", main = "NW2")
for (i in 2:length(nw_features)){
	lines(x = ps, y = nw2_mat[i,], col = colors[i])
	}
legend("bottomright", legend = nw_features, lty = c(1,1,1,1,1,1), col= colors)

plot(x = ps, y = nw3_mat[1,], col = colors[1], type = "l", ylim = c(0,1), xlab = "Sample Size", ylab = "Mean Kendall Rank Correlation", main = "NW3")
for (i in 2:length(nw_features)){
	lines(x = ps, y = nw3_mat[i,], col = colors[i])
	}
legend("bottomright", legend = nw_features, lty = c(1,1,1,1,1,1), col= colors)

plot(x = ps, y = nw4_mat[1,], col = colors[1], type = "l", ylim = c(0,1), xlab = "Sample Size", ylab = "Mean Kendall Rank Correlation", main = "NW4")
for (i in 2:length(nw_features)){
	lines(x = ps, y = nw4_mat[i,], col = colors[i])
	}
legend("bottomright", legend = nw_features, lty = c(1,1,1,1,1,1), col= colors)

plot(x = ps, y = nw5_mat[1,], col = colors[1], type = "l", ylim = c(0,1), xlab = "Sample Size", ylab = "Mean Kendall Rank Correlation", main = "NW5")
for (i in 2:length(nw_features)){
	lines(x = ps, y = nw5_mat[i,], col = colors[i])
	}
legend("bottomright", legend = nw_features, lty = c(1,1,1,1,1,1), col= colors)

plot(x = ps, y = nw6_mat[1,], col = colors[1], type = "l", ylim = c(0,1), xlab = "Sample Size", ylab = "Mean Kendall Rank Correlation", main = "NW6")
for (i in 2:length(nw_features)){
	lines(x = ps, y = nw6_mat[i,], col = colors[i])
	}
legend("bottomright", legend = nw_features, lty = c(1,1,1,1,1,1), col= colors)

plot(x = ps, y = nw7_mat[1,], col = colors[1], type = "l", ylim = c(0,1), xlab = "Sample Size", ylab = "Mean Kendall Rank Correlation", main = "NW7")
for (i in 2:length(nw_features)){
	lines(x = ps, y = nw7_mat[i,], col = colors[i])
	}
legend("bottomright", legend = nw_features, lty = c(1,1,1,1,1,1), col= colors)

plot(x = ps, y = nw8_mat[1,], col = colors[1], type = "l", ylim = c(0,1), xlab = "Sample Size", ylab = "Mean Kendall Rank Correlation", main = "NW8")
for (i in 2:length(nw_features)){
	lines(x = ps, y = nw8_mat[i,], col = colors[i])
	}
legend("bottomright", legend = nw_features, lty = c(1,1,1,1,1,1), col= colors)

## And with varying line type and starker colors:


#### SI Figure 9:

colors <- c("hotpink", "purple", "turquoise", "blue", "forestgreen", "sienna1")

pdf("SIf9a.pdf")
par(mfrow = c(3,3))
plot(x = ps, y = nw1_mat[1,], col = colors[1], type = "l", ylim = c(0,1), xlab = "Sample Size", ylab = "Mean Kendall Rank Correlation", main = "NW1")
for (i in 2:length(nw_features)){
	lines(x = ps, y = nw1_mat[i,], col = colors[i], lty=i)
	}
legend("bottomright", legend = nw_features, lty = c(1:6), col= colors)
dev.off()

pdf("SIf9b.pdf")	
plot(x = ps, y = nw2_mat[1,], col = colors[1], type = "l", ylim = c(0,1), xlab = "Sample Size", ylab = "Mean Kendall Rank Correlation", main = "NW2")
for (i in 2:length(nw_features)){
	lines(x = ps, y = nw2_mat[i,], col = colors[i], lty=i)
	}
legend("bottomright", legend = nw_features, lty = c(1:6), col= colors)
dev.off()


pdf("SIf9c.pdf")
plot(x = ps, y = nw3_mat[1,], col = colors[1], type = "l", ylim = c(0,1), xlab = "Sample Size", ylab = "Mean Kendall Rank Correlation", main = "NW3")
for (i in 2:length(nw_features)){
	lines(x = ps, y = nw3_mat[i,], col = colors[i], lty=i)
	}
legend("bottomright", legend = nw_features, lty = c(1:6), col= colors)
dev.off()

pdf("SIf9d.pdf")
plot(x = ps, y = nw4_mat[1,], col = colors[1], type = "l", ylim = c(0,1), xlab = "Sample Size", ylab = "Mean Kendall Rank Correlation", main = "NW4")
for (i in 2:length(nw_features)){
	lines(x = ps, y = nw4_mat[i,], col = colors[i], lty=i)
	}
legend("bottomright", legend = nw_features, lty = c(1:6), col= colors)
dev.off()

pdf("SIf9e.pdf")
plot(x = ps, y = nw5_mat[1,], col = colors[1], type = "l", ylim = c(0,1), xlab = "Sample Size", ylab = "Mean Kendall Rank Correlation", main = "NW5")
for (i in 2:length(nw_features)){
	lines(x = ps, y = nw5_mat[i,], col = colors[i], lty=i)
	}
legend("bottomright", legend = nw_features, lty = c(1:6), col= colors)
dev.off()

pdf("SIf9f.pdf")
plot(x = ps, y = nw6_mat[1,], col = colors[1], type = "l", ylim = c(0,1), xlab = "Sample Size", ylab = "Mean Kendall Rank Correlation", main = "NW6")
for (i in 2:length(nw_features)){
	lines(x = ps, y = nw6_mat[i,], col = colors[i], lty=i)
	}
legend("bottomright", legend = nw_features, lty = c(1:6), col= colors)
dev.off()

pdf("SIf9g.pdf")
plot(x = ps, y = nw7_mat[1,], col = colors[1], type = "l", ylim = c(0,1), xlab = "Sample Size", ylab = "Mean Kendall Rank Correlation", main = "NW7")
for (i in 2:length(nw_features)){
	lines(x = ps, y = nw7_mat[i,], col = colors[i], lty=i)
	}
legend("bottomright", legend = nw_features, lty = c(1:6), col= colors)
dev.off()

pdf("SIf9h.pdf")
plot(x = ps, y = nw8_mat[1,], col = colors[1], type = "l", ylim = c(0,1), xlab = "Sample Size", ylab = "Mean Kendall Rank Correlation", main = "NW8")
for (i in 2:length(nw_features)){
	lines(x = ps, y = nw8_mat[i,], col = colors[i], lty=i)
	}
legend("bottomright", legend = nw_features, lty = c(1:6), col= colors)
dev.off()

## With labels for article text:
## Article Figure 1:

colors <- c("black", "black", "black", "black", "black", "black")

pdf("f1.pdf")
par(mfrow=c(2,2))

plot(x = ps, y = nw1_mat[1,], col = colors[1], type = "l", ylim = c(0,1), xlab = "Sample Size", ylab = "Mean Kendall Rank Correlation", main = "Network A, Undirected")
for (i in 2:length(nw_features)){
	lines(x = ps, y = nw1_mat[i,], col = colors[i], lty=i)
	}
legend("bottomright", legend = nw_features, lty = c(1:6), col= colors)


plot(x = ps, y = nw2_mat[1,], col = colors[1], type = "l", ylim = c(0,1), xlab = "Sample Size", ylab = "Mean Kendall Rank Correlation", main = "Network B, Undirected")
for (i in 2:length(nw_features)){
	lines(x = ps, y = nw2_mat[i,], col = colors[i], lty=i)
	}
legend("bottomright", legend = nw_features, lty = c(1:6), col= colors)


plot(x = ps, y = nw7_mat[1,], col = colors[1], type = "l", ylim = c(0,1), xlab = "Sample Size", ylab = "Mean Kendall Rank Correlation", main = "Network C, Directed")
for (i in 2:length(nw_features)){
	lines(x = ps, y = nw7_mat[i,], col = colors[i], lty=i)
	}
legend("bottomright", legend = nw_features, lty = c(1:6), col= colors)

plot(x = ps, y = nw8_mat[1,], col = colors[1], type = "l", ylim = c(0,1), xlab = "Sample Size", ylab = "Mean Kendall Rank Correlation", main = "Network D, Directed")
for (i in 2:length(nw_features)){
	lines(x = ps, y = nw8_mat[i,], col = colors[i], lty=i)
	}
legend("bottomright", legend = nw_features, lty = c(1:6), col= colors)
dev.off()


